Using the idea of coverage that was introduced with simulation studies, we can assess the coverage of bootstrapped confidence intervals.

Consider a random sample from a standard normal distribution.

x

 [1] -0.169417734 -0.618778065  0.889289947  0.743545975  0.258047445 -0.007841503
 [7] -1.304085671 -1.842335987 -0.277600642 -1.612780873  0.511142431  0.628909857
[13] -0.361393835 -0.214121337  0.772547387  0.433494812  0.023510859  0.168546163
[19] -0.161329600  0.264790641  0.493521998 -0.307056166 -1.365778025  0.426529696
[25] -0.562316866  1.654235108 -1.086048084 -0.227679632 -1.468276619  0.783348415
[31]  0.190550082  0.519711800 -0.136788620  0.500098666  0.102446975 -0.763293537
[37]  0.334649480 -0.853774135 -0.862946417 -1.064681908  0.583762658 -0.049363852
[43] -0.190355104 -0.862702244 -0.603235304  0.121753100 -0.045877444 -0.029367173
[49]  0.372187947 -0.941077258

Given this, we want to estimate the 25th quantile value.

quantile(x, 0.25)

       25% 
-0.6148924

We can build a confidence interval around this estimate using the bootstrap method.

quantile(bs_samples, c(0.025, 0.975))

      2.5%      97.5% 
-0.9410773 -0.2139728 
bootstrap_ci(x, function(x) quantile(x, 0.25))

      2.5%      97.5% 
-0.9410773 -0.2182428

We also can calculate the actual value of the 25th quantile of the standard normal distribution.

(truth <- qnorm(0.25))

[1] -0.6744898

Given that we know the actual value, we can use a simulation study to assess the coverage of our bootstrap confidence interval.

rr library(tictoc) n_reps <- 1000

tictoc::tic()
bs_ci_values <- replicate(n_reps, bootstrap_ci(rnorm(n), function(x) quantile(x, 0.25)))
tictoc::toc()
83.222 sec elapsed
bs_ci_values[,1:5]

             [,1]       [,2]       [,3]       [,4]       [,5]
2.5%  -0.82931940 -0.9722685 -1.2333750 -0.8899601 -1.3250445
97.5% -0.05235155 -0.5857609 -0.3119804 -0.4367608 -0.2188246

Notice that the above code takes a long time to run. Why?

We can speed things up by running our code in parallel! First, we need to install a couple of useful packages, future.apply and furrr. future.apply provides parallelized versions of the R apply functions while furrr provides parallelized versions of the purrr map functions.

install.packages(c("future.apply", "furrr"))
library(future.apply)
Loading required package: future
library(furrr)

Now, we need to setup our computer to take advantage of the cores we have available.

availableCores()

system 
     8

Run the same code as above, but this time in parallel using future_replicate().

tic()
bs_ci_values <- future_replicate(n_reps, bootstrap_ci(rnorm(n), function(x) quantile(x, 0.25)))
toc()
22.613 sec elapsed
dim(bs_ci_values)

[1]    2 1000
tic()
bs_ci_list <- future_map(1:n_reps, 
                         function(x) bootstrap_ci(rnorm(n), function(x) quantile(x, 0.25)))
toc()
23.416 sec elapsed
length(bs_ci_list)

[1] 1000
head(bs_ci_list)

[[1]]
      2.5%      97.5% 
-0.8573536 -0.1247388 

[[2]]
     2.5%     97.5% 
-1.351951 -0.439162 

[[3]]
      2.5%      97.5% 
-0.8703554 -0.2306619 

[[4]]
       2.5%       97.5% 
-1.01540445 -0.03089471 

[[5]]
      2.5%      97.5% 
-0.8405373 -0.4936884 

[[6]]
      2.5%      97.5% 
-1.1808770 -0.6686824

Now that we have these confidence interval values, we can estimate coverage by determining what proportion of intervals contain the true value.

(coverage <- mean(purrr::map_lgl(bs_ci_list, ~.[1] <= truth & .[2] >= truth)))

[1] 0.939
(coverage <- mean(purrr::map_lgl(bs_ci_list, ~prod(. - truth) < 0)))

[1] 0.939
coverage_ci <- coverage + c(-1, 1)*qnorm(0.975)*sqrt(coverage*(1 - coverage)/n_reps)
c(lower = coverage_ci[1],
  coverage = coverage,
  upper = coverage_ci[2])
    lower  coverage     upper 
0.9241664 0.9390000 0.9538336
LS0tCnRpdGxlOiAiQm9vdHN0cmFwIENvdmVyYWdlIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpgYGB7ciBzZXR1cCwgaW5jbHVkZSA9IEZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQoY2FjaGUgPSBUUlVFLCBtZXNzYWdlID0gRkFMU0UpCmBgYAoKVXNpbmcgdGhlIGlkZWEgb2YgY292ZXJhZ2UgdGhhdCB3YXMgaW50cm9kdWNlZCB3aXRoIHNpbXVsYXRpb24gc3R1ZGllcywgd2UgY2FuIAphc3Nlc3MgdGhlIGNvdmVyYWdlIG9mIGJvb3RzdHJhcHBlZCBjb25maWRlbmNlIGludGVydmFscy4KCkNvbnNpZGVyIGEgcmFuZG9tIHNhbXBsZSBmcm9tIGEgc3RhbmRhcmQgbm9ybWFsIGRpc3RyaWJ1dGlvbi4KCmBgYHtyfQpuIDwtIDUwCnggPC0gcm5vcm0obikKYGBgCgpHaXZlbiB0aGlzLCB3ZSB3YW50IHRvIGVzdGltYXRlIHRoZSAyNXRoIHF1YW50aWxlIHZhbHVlLgpgYGB7cn0KcXVhbnRpbGUoeCwgMC4yNSkKYGBgCgpXZSBjYW4gYnVpbGQgYSBjb25maWRlbmNlIGludGVydmFsIGFyb3VuZCB0aGlzIGVzdGltYXRlIHVzaW5nIHRoZSBib290c3RyYXAgbWV0aG9kLgpgYGB7cn0KYnNfc2FtcGxlcyA8LSByZXBsaWNhdGUoMTAwMCwgcXVhbnRpbGUoc2FtcGxlKHgsIHJlcGxhY2UgPSBUUlVFKSwgMC4yNSkpCmhpc3QoYnNfc2FtcGxlcykKcXVhbnRpbGUoYnNfc2FtcGxlcywgYygwLjAyNSwgMC45NzUpKQoKYm9vdHN0cmFwX2NpIDwtIGZ1bmN0aW9uKHgsIGYsIGFscGhhID0gMC4wNSwgbl9zYW1wbGVzID0gMTAwMCkgewogIGJzX291dCA8LSByZXBsaWNhdGUobl9zYW1wbGVzLCBmKHNhbXBsZSh4LCByZXBsYWNlID0gVFJVRSkpKQogIHF1YW50aWxlKGJzX291dCwgYyhhbHBoYS8yLCAxIC0gYWxwaGEvMikpCn0KYGBgCgpgYGB7cn0KYm9vdHN0cmFwX2NpKHgsIGZ1bmN0aW9uKHgpIHF1YW50aWxlKHgsIDAuMjUpKQpgYGAKCldlIGFsc28gY2FuIGNhbGN1bGF0ZSB0aGUgYWN0dWFsIHZhbHVlIG9mIHRoZSAyNXRoIHF1YW50aWxlIG9mIHRoZSBzdGFuZGFyZCBub3JtYWwKZGlzdHJpYnV0aW9uLgpgYGB7cn0KKHRydXRoIDwtIHFub3JtKDAuMjUpKQpgYGAKCkdpdmVuIHRoYXQgd2Uga25vdyB0aGUgYWN0dWFsIHZhbHVlLCB3ZSBjYW4gdXNlIGEgc2ltdWxhdGlvbiBzdHVkeSB0byBhc3Nlc3MgdGhlCmNvdmVyYWdlIG9mIG91ciBib290c3RyYXAgY29uZmlkZW5jZSBpbnRlcnZhbC4KCmBgYHtyfQpsaWJyYXJ5KHRpY3RvYykKbl9yZXBzIDwtIDEwMDAKYGBgCgpgYGB7cn0KdGljdG9jOjp0aWMoKQpic19jaV92YWx1ZXMgPC0gcmVwbGljYXRlKG5fcmVwcywgYm9vdHN0cmFwX2NpKHJub3JtKG4pLCBmdW5jdGlvbih4KSBxdWFudGlsZSh4LCAwLjI1KSkpCnRpY3RvYzo6dG9jKCkKYGBgCgpgYGB7cn0KZGltKGJzX2NpX3ZhbHVlcykKYnNfY2lfdmFsdWVzWywxOjVdCmBgYAoKTm90aWNlIHRoYXQgdGhlIGFib3ZlIGNvZGUgdGFrZXMgYSAqbG9uZyogdGltZSB0byBydW4uIFdoeT8KCldlIGNhbiBzcGVlZCB0aGluZ3MgdXAgYnkgcnVubmluZyBvdXIgY29kZSBpbiBwYXJhbGxlbCEgRmlyc3QsIHdlIG5lZWQgdG8gaW5zdGFsbAphIGNvdXBsZSBvZiB1c2VmdWwgcGFja2FnZXMsIGBmdXR1cmUuYXBwbHlgIGFuZCBgZnVycnJgLiBgZnV0dXJlLmFwcGx5YCBwcm92aWRlcwpwYXJhbGxlbGl6ZWQgdmVyc2lvbnMgb2YgdGhlIFIgYGFwcGx5YCBmdW5jdGlvbnMgd2hpbGUgYGZ1cnJyYCBwcm92aWRlcyBwYXJhbGxlbGl6ZWQKdmVyc2lvbnMgb2YgdGhlIGBwdXJycmAgbWFwIGZ1bmN0aW9ucy4KCmBgYHtyLCBldmFsPUZBTFNFfQppbnN0YWxsLnBhY2thZ2VzKGMoImZ1dHVyZS5hcHBseSIsICJmdXJyciIpKQpgYGAKCmBgYHtyfQpsaWJyYXJ5KGZ1dHVyZS5hcHBseSkKbGlicmFyeShmdXJycikKYGBgCgpOb3csIHdlIG5lZWQgdG8gc2V0dXAgb3VyIGNvbXB1dGVyIHRvIHRha2UgYWR2YW50YWdlIG9mIHRoZSBjb3JlcyB3ZSBoYXZlIGF2YWlsYWJsZS4KCmBgYHtyfQphdmFpbGFibGVDb3JlcygpCnBsYW4oIm11bHRpcHJvY2VzcyIpCmBgYAoKIyBSdW4gdGhlIHNhbWUgY29kZSBhcyBhYm92ZSwgYnV0IHRoaXMgdGltZSBpbiBwYXJhbGxlbCB1c2luZyBgZnV0dXJlX3JlcGxpY2F0ZSgpYC4KCmBgYHtyfQp0aWMoKQpic19jaV92YWx1ZXMgPC0gZnV0dXJlX3JlcGxpY2F0ZShuX3JlcHMsIGJvb3RzdHJhcF9jaShybm9ybShuKSwgZnVuY3Rpb24oeCkgcXVhbnRpbGUoeCwgMC4yNSkpKQp0b2MoKQpgYGAKCmBgYHtyfQpkaW0oYnNfY2lfdmFsdWVzKQpgYGAKCgpgYGB7cn0KdGljKCkKYnNfY2lfbGlzdCA8LSBmdXR1cmVfbWFwKDE6bl9yZXBzLCAKICAgICAgICAgICAgICAgICAgICAgICAgIGZ1bmN0aW9uKHgpIGJvb3RzdHJhcF9jaShybm9ybShuKSwgZnVuY3Rpb24oeCkgcXVhbnRpbGUoeCwgMC4yNSkpKQp0b2MoKQpgYGAKCmBgYHtyfQpsZW5ndGgoYnNfY2lfbGlzdCkKYGBgCgpgYGB7cn0KaGVhZChic19jaV9saXN0KQpgYGAKCgpOb3cgdGhhdCB3ZSBoYXZlIHRoZXNlIGNvbmZpZGVuY2UgaW50ZXJ2YWwgdmFsdWVzLCB3ZSBjYW4gZXN0aW1hdGUgY292ZXJhZ2UgYnkKZGV0ZXJtaW5pbmcgd2hhdCBwcm9wb3J0aW9uIG9mIGludGVydmFscyBjb250YWluIHRoZSB0cnVlIHZhbHVlLgoKYGBge3J9Cihjb3ZlcmFnZSA8LSBtZWFuKHB1cnJyOjptYXBfbGdsKGJzX2NpX2xpc3QsIH4uWzFdIDw9IHRydXRoICYgLlsyXSA+PSB0cnV0aCkpKQpgYGAKCmBgYHtyfQooY292ZXJhZ2UgPC0gbWVhbihwdXJycjo6bWFwX2xnbChic19jaV9saXN0LCB+cHJvZCguIC0gdHJ1dGgpIDwgMCkpKQpgYGAKCgpgYGB7cn0KY292ZXJhZ2VfY2kgPC0gY292ZXJhZ2UgKyBjKC0xLCAxKSpxbm9ybSgwLjk3NSkqc3FydChjb3ZlcmFnZSooMSAtIGNvdmVyYWdlKS9uX3JlcHMpCmMobG93ZXIgPSBjb3ZlcmFnZV9jaVsxXSwKICBjb3ZlcmFnZSA9IGNvdmVyYWdlLAogIHVwcGVyID0gY292ZXJhZ2VfY2lbMl0pCmBgYAo=